library(foreign)
library(DAMisc)


cand <- read.csv('prq_fig4_5.csv')



#############################################################################
#FIGURE 4: Comparison of Model Predictions, All Cases and Marginal Districts#
#############################################################################

N <- max(cand$decade, na.rm = TRUE)

###############
#pooled models#
###############

p1 <- array(NA, dim = c(N,1))
p2 <- array(NA, dim = c(N,1))
p3 <- array(NA, dim = c(N,1))


for(i in 1:N){
  m1 <- glm(dwin ~   dpres , subset(cand, decade == i ), family = binomial)
  p1[i] <- unlist(pre(m1)$pcp)
  
  m2 <- glm(dwin ~   dpres + inc_t, subset(cand, decade == i ), family = binomial)
  p2[i] <- unlist(pre(m2)$pcp)
  
  #m3 <- glm(dwin ~  party, subset(cand, yearID == i), family = binomial)
  #p3[i] <- p2[i] - p1[i] 
  
}


p3 <- unlist(p2) - unlist(p1)
#Plot comparison of PRE


dec.labs <- c("1840s", "1850s", "1860s", "1870s", "1880s", "1890s", "1900s", 
              "1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", 
              "1990s", "2000s", "2014-18")
prop.labs <- c('0.0', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', 
               '0.7', '0.8', '0.9', '1.0')
dec <- 1:N

#two figs per plt
par(mfrow = c(1,2))

#plot G-K over time
par(mar = c(5, 4, 2, 2))

plot(dec, p1, ylim = c(0,1), type = 'n', xlab='', ylab='', xaxt='n', yaxt='n' )


#code puts axis text and an angle
axis(side = 1, at = dec, labels = FALSE ) 
text(dec, par("usr")[3] - 0.025, labels = dec.labs, srt = 45, pos = 1, xpd = TRUE, cex = 0.8)

axis(side = 2, at = seq(0, 1, 0.2),  las = 2)
#segments(x0 = tick, x1 = tick, y0 = inc.lo, y1 = inc.hi, lwd =2)
lines(dec, p1, col = "grey", lwd = 3)
lines(dec, p3, lwd = 3,  lty = 2)
#lines(tick, p3, lwd = 3,  lty = 1)


mtext("Year", side = 1, line = 2.75)
mtext("Percent Correctly Predicted", side = 2, line = 2.75)
mtext("All Distircts", side = 3, line = 0.5, cex = 1.25)

legend('topleft', bty='n', legend = c('Model 1', 'Model 2 - Model 1'), lwd = c(3,3), col = c('grey', 'black'), cex = 0.8, lty=c(1,2))




################################
#marginal pres districts models#
################################

p1 <- array(NA, dim = c(N,1))
p2 <- array(NA, dim = c(N,1))
p3 <- array(NA, dim = c(N,1))


for(i in 1:N){
  m1 <- glm(dwin ~   dpres , subset(cand, dpres > 45 & dpres < 55 & decade == i ), family = binomial)
  p1[i] <- unlist(pre(m1)$pcp)
  
  m2 <- glm(dwin ~   dpres + inc_t, subset(cand, dpres > 45 & dpres < 55 & decade == i ), family = binomial)
  p2[i] <- unlist(pre(m2)$pcp)
  
  #m3 <- glm(dwin ~  party, subset(cand, yearID == i), family = binomial)
  #p3[i] <- p2[i] - p1[i] 
  
}


p3 <- unlist(p2) - unlist(p1)
#Plot comparison of PRE


dec.labs <- c("1840s", "1850s", "1860s", "1870s", "1880s", "1890s", "1900s", 
              "1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", 
              "1990s", "2000s", "2014-18")
prop.labs <- c('0.0', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', 
               '0.7', '0.8', '0.9', '1.0')
dec <- 1:N

#plot G-K over time
par(mar = c(5, 4, 2, 2))

plot(dec, p1, ylim = c(0,1), type = 'n', xlab='', ylab='', xaxt='n', yaxt='n' )


#code puts axis text and an angle
axis(side = 1, at = dec, labels = FALSE ) 
text(dec, par("usr")[3] - 0.025, labels = dec.labs, srt = 45, pos = 1, xpd = TRUE, cex = 0.8)

axis(side = 2, at = seq(0, 1, 0.2),  las = 2)
#segments(x0 = tick, x1 = tick, y0 = inc.lo, y1 = inc.hi, lwd =2)
lines(dec, p1, col = "grey", lwd = 3)
lines(dec, p3, lwd = 3,  lty = 2)
#lines(tick, p3, lwd = 3,  lty = 1)


mtext("Decade", side = 1, line = 2.75)
mtext("Percent Correctly Predicted", side = 2, line = 2.75)
mtext("Marginal Districts", side = 3, line = 0.5, cex = 1.25)

legend('topleft', bty='n', legend = c('Model 1', 'Model 2 - Model 1'), lwd = c(3,3), col = c('grey', 'black'), cex = 0.8, lty=c(1,2))


###################################################################################
#FIGURE 5:Comparison of Model Predictions, Presidential and Midterm Election Years#
###################################################################################

par(mfrow = c(1,1))

#################
#midterm  models#
#################

p1m <- array(NA, dim = c(N,1))


for(i in 1:N){
  m1 <- glm(dwin ~   dpres , subset(cand, pres_elec == 0 & decade == i ), family = binomial)
  p1m[i] <- unlist(pre(m1)$pcp)
  
  
}


###################
#pres year  models#
###################

p1p <- array(NA, dim = c(N,1))

for(i in 1:N){
  m1 <- glm(dwin ~   dpres , subset(cand, pres_elec == 1 & decade == i ), family = binomial)
  p1p[i] <- unlist(pre(m1)$pcp)
}



########################################
#COMPARISON OF MIDTERM AND PRESIDENTIAL#
########################################

dec.labs <- c("1840s", "1850s", "1860s", "1870s", "1880s", "1890s", "1900s", 
              "1910s", "1920s", "1930s", "1940s", "1950s", "1960s", "1970s", "1980s", 
              "1990s", "2000s", "2014-18")
prop.labs <- c('0.0', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', 
               '0.7', '0.8', '0.9', '1.0')
dec <- 1:N

#plot G-K over time
par(mar = c(5, 4, 2, 2))

plot(dec, p1, ylim = c(0,1), type = 'n', xlab='', ylab='', xaxt='n', yaxt='n' )


#code puts axis text and an angle
axis(side = 1, at = dec, labels = FALSE ) 
text(dec, par("usr")[3] - 0.025, labels = dec.labs, srt = 45, pos = 1, xpd = TRUE, cex = 0.8)

axis(side = 2, at = seq(0, 1, 0.2),  las = 2)
#segments(x0 = tick, x1 = tick, y0 = inc.lo, y1 = inc.hi, lwd =2)
lines(dec, p1p, col = "grey", lwd = 3)
lines(dec, p1m, lwd = 3,  lty = 2)
#lines(tick, p3, lwd = 3,  lty = 1)


mtext("Decade", side = 1, line = 2.75)
mtext("Percent Correctly Predicted", side = 2, line = 2.75)
#mtext("Presidential", side = 3, line = 0.5, cex = 1.25)

legend('bottomleft', bty='n', legend = c('Presidential', 'Midterm'), lwd = c(3,3), col = c('grey', 'black'), cex = 0.8, lty=c(1,2))
